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Self-organizing nanocheckerboards have been experimentally fabricated in Mn-based spinels, but 
have not yet been explained with first principles. Using density-functional-theory, we explain the 
phase diagram of the ZnMnxGa 2 -x 04 system and the origin of nanocheckerboards. We predict total 
phase separation at zero temperature, then show the combination of kinetics, thermodynamics, and 
Jahn-Teller physics that generates the system’s observed behavior. We find the {011} surfaces 
are strongly-preferred energetically, which mandates checkerboard ordering by purely geometrical 
considerations. 


Experimental observation demonstrates intriguing 
nanoscale compositional ordering in a variety of mate¬ 
rial alloys. These include noble-metal-alloy nanochecker¬ 
boards m, BaTi03-CoFe2 04 nanopillars [1], and an 
assortment of manganite-spinel nanocheckerboards 0- 
[9]. Nanoscale phenomena are inherently difficult to treat 
with quantum mechanics’ first principles, due to the pro¬ 
hibitive scaling of electronic-structure methods. Previ¬ 
ous theoretical studies [Mm used phase-field models 
[m nn to simulate nanocheckerboard formation. How¬ 
ever, those models rely upon coefficients chosen without 
first-principles justification. In contrast, our work reveals 
the origin of nanocheckerboards from first-principles. 

Here we examine the experimentally well-characterized 
manganite spinels These cooperative- 

Jahn-Teller crystals, upon doping with certain transi¬ 
tion metals, organize into nanocheckerboards. Exper¬ 
iments showed that high-temperature mixing, followed 
by slow cooling, yields a spontaneously-formed checker¬ 
board whose squares alternate tetragonal Mn-rich and 
cubic Mn-poor phases. These checkerboards emerge from 
the cross-section of spontaneously-aligned nanorods. Yeo 
et al. |5] fabricated self-assembling nanocheckerboards 
from ZnGa204 (ZGO) -I- ZnMn204 (ZMO), comprised 
of ^ 4 X 4 X 70 nm^ nanorods. Later work grew 
checkerboards with nanorods over 700 nm long on a 
MgO substrate n [7]. Gheckerboards were later ex¬ 
tended to other Mn-based spinels, first MgMnxFe2-x04 
(MMFO) |8] and then the tunable-sized checkerboards of 
Goo,6MnxFe2,4-x04 (GMFO) 0. The latter are notable 
for a patterning that alternates ferro- and paramagnetic 
phases, which yields potential for ultrahigh-density infor¬ 
mation storage. 

We analyze the nanocheckerboard system 
ZnMnxGa2-x04, whose relative simplicity renders 
it a minimal prototype. We sketch the main experimen¬ 
tal results of Ref. |S]in Figure ZGO is a cubic spinel 
while ZMO is a tetragonal spinel (c/a = 1.14 [H]): 
which immediately leads to some anomalous differences 
between experiment and naive expectations. First, 
room-temperature experiments reveal solid solutions 
at non-negligible concentrations, despite the disparate 
crystal structures of the end-members, in violation 



FIG. 1. Sketch of experimental results presented in Ref. O 
For a;Mn < 0.25 a single-phase cubic structure (G) appears, 
while for a^Mn > 0.85 it is a single-phase tetragonal structure 
(T). For intermediate concentrations, a single-phase tetrago¬ 
nal structure (T') is observed upon rapid quenching (the high- 
temperature phase), while slower cooling generates phase- 
separated nanocheckerboards with (Oil) and (011) inter¬ 
faces. Arrows in the tetragonal regions show the tetragonally- 
elongated direction. The concentrations on the abscissa are 
those for which data were reported in Ref. O The concen¬ 
trations within the nanocheckerboards were not reported as 
being measured directly. 


of the Hume-Rothery rules. Second, x-ray diffraction 
experiments show a cubic structure up to 25% ZMO 
(and 75% ZGO), whereas a non-negligible tetragonality 
would be expected at this concentration, regardless 
of the origin of the observed solubility. Third, in the 
region where ZMGO phase-separates, experiment shows 
checkerboard formations instead of traditional spinodal 
decomposition. 

To address these issues, we compute the energetics of 
ZnMnxGa2_x04 using density functional theory as im¬ 
plemented in the Vienna Ab-initio Simulation Package 
(VASP) [TSHII]. We use a plane-wave cutoff energy 
of 415 eV, the PW91 generalized gradient approxima¬ 
tion (GGA) functional [221123], and projector augmented 
wave (PAW) based pseudopotentials [23]. All calcula¬ 
tions are initialized with ferromagnetic ordering for sim¬ 
plicity, a well-justihed approach for the Neel tempera¬ 
ture of merely « 20 K [mnauB]. We compute a phase 
diagram via a cluster expansion [27| of the DFT ener¬ 
getics with Monte-Carlo simulations, as implemented in 
the Alloy Theoretic Automated Toolkit (ATAT) package 
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Both ZGO and ZMO crystallize in the spinel struc¬ 
ture, with nominal valences of Zn^+(Ga, Mn) 2 ~'' 04 “. Zn 
occupies the tetrahedral (“A”) sites and Ga / Mn oc¬ 
cupy the octahedral (“B”) sites, with negligible inversion 
[351 [33]. Therefore, our study focuses on the effect of Ga 
/ Mn occupation of the B sites. Each B-centered octahe¬ 
dron shares edges with six others, which couples their an¬ 
ionic distortions. ZGO forms a cubic spinel (space group 
FdZm) [33], while ZMO is a tetragonal spinel (space 
group lAi/amd) with significant distortion c/a = 1.14 
m- The T' of Figure (high-temperature fully-mixed 
ZMGO) has c/a « 1.06 [5]. We have provided lattice 
constants from the literature, and their DFT-calculated 
analogs, in Supplementary Material (Table SI); DFT cal¬ 
culations agree with experiment. ZMO’s sizable tetrag¬ 
onal distortion is due to the Jahn-Teller (JT) effect in 
the Mn^'*’ ions, where the configuration in a high-spin 
octahedral environment causes the Cg orbitals to break 
symmetry by a tetragonal distortion. This orbital order¬ 
ing leads to a martensitic cubic-to-tetragonal transition 
in a variety of crystals, including spinels [31H35] . For con¬ 
sistency, we take [001] to be the JT-distorted direction. 

Experiments, summarized in Figure[^ show solid solu¬ 
tions for XMn < 0.25 and Xmu > 0.85. Yet the cubic and 
strongly-tetragonal crystal structures of the end mem¬ 
bers are expected to be immiscible: Placing a non-JT 
octahedron in a tetragonal JT-active environment, and 
vice versa, costs energy. Our DFT calculations quanti¬ 
tatively verify the qualitative Hume-Rothery rules: We 
generated 192 supercells (< 42 atoms) then fully relaxed 
the structures. (See Supplementary Material for further 
calculation details.) Figure [^a) shows their formation 
energies, demonstrating that mixing the disparate crys¬ 
tal structures incurs an energy cost [39]. A ground- 
state search using the cluster expansion (via the afore¬ 
mentioned ATAT package) revealed no lower-energy su¬ 
percells. Therefore, the zero-temperature ground state 
of ZMGO is actually phase separation into bulk ZGO 
and ZMO, a conclusion absent from previous non-first- 
principles analyses [IIHIll. There is no chemical or phys¬ 
ical reason to believe that ZGO and ZMO should mix at 
anything but elevated temperatures. 

To further verify this, we computed a phase diagram 
via a cluster expansion, using the aforementioned ATAT 
package. (See Supplementary Material for extensive cal¬ 
culation details.) As shown in Figure [^b), the zero- 
temperature stable phases are immiscible bulk ZGO and 
ZMO. In fact, the observed miscibility at XMn £ 0.25 and 
a:Mn > 0.85 is unreasonable at all but extreme tempera¬ 
tures. 

We therefore attribute the anomalous miscibility to ki¬ 
netic limitations, i.e. diffusion constraints. This is cor¬ 
roborated by the experimental observation that slower 
cooling (i.e. better diffusion) leads to larger checkers 
(i.e. less miscibility) [51[n]. In fact, kinetics are the ob- 
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FIG. 2. (a) Energies of formation of 192 ZnMnxGa 2 -x 04 

supercells, (b) Phase diagram of mixed ZGO and ZMO, us¬ 
ing ATAT’s EMG2 code. The stable phases (with respect 
to a semi-grand-canonical ensemble, where a;Mn may vary) 
are shown at varying temperatures. The phases are marked: 
Gubic (G), tetragonal (T), phase-separated (C-l-T), and high- 
temperature fully-mixed tetragonal (T'). 


vious origin of the miscible T' state observed upon rapid 
quenching (shown in Figure [^. DFT’s prediction for the 
T' state for xun = 0.5 (taken as the mean of all calcu¬ 
lated structures) agrees well with the experimental mea¬ 
surements. For example, DFT predicts c/a = 1.07, com¬ 
parable with experiment’s c/a k, 1.06. (See Supplemen¬ 
tary Material for all data.) Similarly, kinetic limitations 
must cause the apparent solubility of XMn < 0.25 and 
a^Mn > 0.85. For example, if diffusion essentially freezes 
by e.g. 900 K, the system cannot separate into 100/0% 
Mn mixtures and will remain a high-entropy frozen solid 
solution. Unfortunately, quantitatively predicting kinet¬ 
ics (including checker size) requires detailed knowledge of 
the diffusion mechanisms, coupled with complex JT lat¬ 
tice dynamics, which is beyond the scope of this Letter. 

Having explained that the apparent solubility is due 
to kinetics, we must address the second anomaly, the 
experimental observation of a cubic structure for a;Mn < 
0.25. In contrast, DFT predicts a tetragonal structure 
(c/a « 1.03 for all calculated structures at XMn = 0.25) 
due to the cooperative Jahn-Teller effect. 

We propose that the cubic structure of the Mn-poor 
phase is caused by noncooperative JT distortions at finite 
temperatures. It is well-known that, as in many spinels, 
ZMO undergoes a phase transition to a cubic spinel above 
approximately 1323 K, due to noncooperative JT distor¬ 
tions (each octahedron distorting in a random direction) 
[ 3511331135 ] . The transition temperature scales with dop- 
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ing: Within mean-field theory it is approximately lin¬ 
ear with doping |40j : experimentally, a variety of spinels 
transition at « 9260 (c/a — 1) K [35], where c/o is the 
tetragonal distortion induced by the cooperative distor¬ 
tion m- As illustrated in Figure both these trends 
imply that at room temperature, structures of a; < 0.25 
are cubic due to this transition, while x > 0.25 are tetrag¬ 
onal. Similar “critical concentrations” of JT-ions have 
been documented in the literature [101 HO] • This is also 
consistent with the conclusion of Noh et al. [H] that 
although XMn = 0.25 has an XRD pattern of a cubic 
spinel, it has a rather large JT splitting of ^ 0.7 eV, 
i.e. a high-spin electronic configuration, which leads to 
JT distortion. (See Supplementary Material for detailed 
analysis of the reliability of GGA’s high-spin prediction.) 
Therefore, low-temperature measurements should reveal 
DFT’s tetragonally-distorted structure for x < 0.25. 
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gies indicate stable interface directions [3S|. Figure Qa) 
shows the energies of formation for various slab thick¬ 
nesses. (Larger thicknesses are inaccessible due to large, 
convergence-challenged supercells.) In agreement with 
experiment, DFT prefers the {011} surfaces. The energy 
relative to the next-preferred surface is about 10 meV/B 
(160 meV per cubic unit cell) at only 1.5 nm, and pre¬ 
sumably larger for the 4 nm nanocheckerboards. 


Slab Formation Energies 
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FIG. 3. Phase diagram for the noncooperative Jahn-Teller 
transition in spinels. DFT calcnlates the distortion (left axis) 
of our ZMGO strnctnres (as in Figure[^a)) at 0 K. The crit¬ 
ical temperature for the tetragonal-to-cubic transition (right 
axis) is obtained based on the fitting Tc « 9260(c/o — 1) K, 
empirically accurate for a variety of spinels |35| . Therefore at 
room temperatnre, XMn fS 0.25 is cnbic. 

Therefore, the cubic structure is due to the high- 
entropy noncooperative distortion. We should note that 
these finite-temperature Jahn-Teller lattice dynamics are 
expected to enhance miscibility, as noted by other ex¬ 
periments |44j . For example, whereas experiment places 
the boundaries at 25/85% Mn, our phase diagram (Fig¬ 
ure [^b)) shows that ZMO is more tolerant of Ga-doping 
than ZGO of Mn-doping. This discrepancy is likely due 
to the lack of noncooperative distortions in the cluster- 
expansion model, although a quantitative treatment is 
beyond the scope of this Letter. However, this alone 
would not suffice to cause appreciable solubility near 
room temperature. 

Now we turn to 0.25 < XMn < 0.85, where slowly- 
cooled samples phase-separate on a diffusion-limited 
scale. We seek to explain how this leads to nanochecker¬ 
boards rather than traditional spinodal decomposition. 
We calculate the energy of joining a slab of ZGO to 
a slab of ZMO along a particular surface; lower ener- 


FIG. 4. (a) Formation energies of ZGO/ZMO slabs layered 

in the given direction, normalized per B-site, plotted against 
the average thickness of the ZGO and ZMO slabs, (b) Same, 
multiplied by the number of B-sites (N) and divided by the 
interface area (A), giving a per-area normalization, (c) Fit¬ 
ting cubic and tetragonal domains with a coherent interface. 
Arrows show the tetragonally elongated direction. The {011} 
surfaces of [001]-distorted tetragonal and cubic crystals form 
a coherent interface within the lattice. Due to the bent an¬ 
gles at the interfaces, this can occur only as checkerboards 
{Sx = 0, with a four-corner geometry). The rotation 0 can be 
found analytically (see text). Dashed lines represent lattice 
constants, which must match for coherent interfaces. 

Due to symmetry, only five directions are calculated di¬ 
rectly. We search formation energies of multilayer slabs 
of thickness t oriented in an arbitrary direction k by per¬ 
forming an expansion in the symmetry-adapted (tetrag¬ 
onal) harmonics, a standard group-theory methodology 
[46l Wf\ . It is apparent from Figure |^b) that the to¬ 
tal formation energy can be approximated by £'(k, t) = 
A{ca{^)+tcv{'k)) , where A is the cross-sectional area and 
Ca, Cy are constants. We perform an expansion of Ca, Cy in 
the four lowest-order harmonics via a least-squares fit for 
the five calculated values of k. (Fitting data appear in 
Supplementary Material.) We thus confirm that (Oil), 
and equivalently (Oil), are the lowest-energy surfaces. 
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Hence physics dictates that when our system co¬ 
herently mixes cubic and tetragonal phases, it forms 
(Oil) and (Oil) interfaces. Now pure geometry dictates 
checkerboard configurations for coherent interfaces. As 
shown in Figure |^c), due to bent angles at the inter¬ 
faces, the only configuration that retains a coherent lat¬ 
tice is the checkerboard formation. It was previously 
observed that alternating cubic phases are rotated by 
a few degrees, while alternating tetragonal phases are 
rotated by 90° 0 11 [H]- The origin is now obvious 
from Figure l^c), with the cubic-phase angle of rotation 
9c = 7r/2 — 2tan“^ a/c. This agrees with measured val¬ 
ues within < 1° (data in Supplementary Material, Table 
S2). This analysis relates our checkerboards to the CoPt 
cubic-tetragonal nanostructures [Diin] and other lattice- 
induced interface rotations [15]. 

Whence does the {011} preference originate? Previous 
work showed the importance of both strain and local ionic 
distortions (i.e. k = 0 phonons) to the Jahn-Teller effect 
[49| . Our slabs’ formation energy consists of strain en¬ 
ergy, due to biaxial lattice-matching, and contact energy, 
due to local ionic distortions and local binding energies. 
Specifically, for two slabs of thickness t joined in direction 
k, we write: 

^formation — A^^i^strain (k, t) 4- AF/contact (k, t) (1) 

where the units now ensure that the latter two E terms 
converge for for t —>■ oo. As shown in Figure [^a), the 
strain energy can be computed directly by relaxing ionic 
coordinates for ZGO and ZMO separately, with strained 
lattice vectors (the “strained bulk” calculation). Then 
the contact energy is simply the difference between the 
formation energy and the strain energy. (Detailed equa¬ 
tions appear in Supplementary Material.) 

These energies are shown in Figure [^b-c). We note 
that {111} and especially {011} have negative contact 
energies, meaning the strained layered-slab heterostruc¬ 
ture is more stable than the strained ZGO and ZMO sep¬ 
arated bulk. Hence the local ionic distortions couple ben¬ 
eficially to the lattice strain, i.e. atomic rearrangements 
partially alleviate the energy penalties of lattice strain. 
This is predominantly concentrated in the breathing and 
tetragonal-distortion local modes {qi and ga in the nota¬ 
tion of Ref. HOD- Remarkably, the contact energy does 
not converge in the accessible thicknesses, due to the long 
range of the Jahn-Teller effect. Rather, even a remark¬ 
able distance from the surface, interoctahedron coupling 
leads to intracellular atomic displacement. Therefore the 
{011} preference originates in a beneficial coupling be¬ 
tween strain and long-range atomic displacements (con¬ 
trary to previous non-first-principles work neglecting the 
latter mM)- 

In conclusion, we have presented the physics of 
nanocheckerboards based on first-principles calculations. 
We established that the thermodynamic ground state is 
complete phase separation. The incomplete separation 
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FIG. 5. (a) Three types of ZMGO structures used to 

decompose slab formation energies. Black arrows represent 
cross-sectional biaxial strain, calculated by fully relaxing the 
strained heterostructure. All calculations use periodic bound¬ 
ary conditions, (b-c) Slab formation-energy decomposition, as 
defined in the text. Motivation for unit choice is described in 
the text; for comparison, note that a cubic unit cell has 16 B- 
sites and a (lOO)-cross-sectional area of approximately 74A . 
Negative contact energy, as in the {011} surface, indicates 
energetic preference for the strained layered-slab heterostruc¬ 
ture over the strained bulk. 


originates in diffusion limitations, leading to nanoscale 
phase domains. We explained the observed cubic crys¬ 
tal structure at XMn = 0.25 based on noncooperative 
Jahn-Teller distortions at room temperature. There¬ 
fore, although ZMGO’s ground state is bulk-incoherent, 
the diffusion-limited structure is bulk-coherent (using the 
terminology of Ref. EH). This bulk coherence leads to 
phase separation of cubic and tetragonal phases along 
{011} surfaces, which we showed from first principles. 
This, in the presence of kinetic constraints, automati¬ 
cally leads to checkerboards, based on pure geometry. 
The preference for {011} surfaces originates in benefi¬ 
cial coupling between local distortions and lattice strain. 
Further quantitative understanding will require robust 
models for the Jahn-Teller effect in doped materials at 
finite temperatures. 
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Science User Eacility supported by the Office of Science 
of the U.S. Department of Energy under Contract No. 
DE-AC02-05CH11231. The authors acknowledge sup¬ 
port from a DARPA Young Eaculty Award, Grant No. 
D13AP00051. 
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Supplemental Materials 

Here we describe various calculation details in the Letter. 

• We provide lattice constants for various crystals according to both experiment and theory. 

• We detail the cluster expansion (CE) and related phase-diagram calculations. We discuss the quality of the CE 
and its challenges for quantitative prediction. 

• We justify GGA’s calculation of a high-spin Mn electronic state. 

• We provide the calculation details for our expansion of surface energies in the symmetry-adapted harmonics. 

• We explain more details of the decomposition of multilayer-slab formation energies into strain and contact 
components. 


LATTICE CONSTANTS 

Here we detail the lattice constants of the bulk ZMO, ZGO, and mixed structure (T'). Data are provided in Tablejl] 
DFT accurately captures the lattice constants of ZGO, ZMO, and the high-temperature fully-mixed ZMGO (T'). We 
model T' as the average of all calculated structures. These calculations agree with experiment 01121133] , within the 
1-2% overestimation characteristic of GGA calculations. We adopt the convention where a cubic crystal has c/a = 1, 
rather than 


Lattice param. 

Exp. DFT 

DFT Err. 

ZGO bulk a 

8.34A 8.46A 

1.5% 

ZMO bulk a 

8.09A 8.15A 

0.8% 

ZMGO hi-T a 

8.2A 8.31A 

1.3% 

ZGO->ZMGO 5a 

-1.6% -1.8% 


ZMO^ZMGO 5a 

1.4% 1.9% 


ZGO bulk (a=) c 

8.34A 8.46A 

1.5% 

ZMO bulk c 

9.24A 9.41A 

1.8% 

ZMGO hi-T c 

8.7A 8.89A 

2.2% 

ZGO->ZMGO 5c 

4.4% 5.1% 


ZMO^ZMGO 5c 

-5.9% -5.5% 


ZMO bulk c/a 

1.14 

1.15 

0.9% 

ZMGO hi-T c/a 

1.06 

1.07 

0.9% 


TABLE I. Lattice constants. Experimental data from Ref. 151 1171 and 1331 


Table [ h] compares theory and experiment for the rotation angle 9 of the cubic nanochecker domains. (See the Letter 
for the equation used.) It is apparent that this simple model successfully calculates the rotation to within < 1°. 


CLUSTER EXPANSION CALCULATIONS 

Here we detail the cluster expansion (CE) and related calculations used to calculate the phase diagram in Figure 
l^b). The CE was based on ZMGO structures. These calculations relied on VASP and the ATAT package; references 
appear in the Letter. 


Ref. 

Substrate 

9 meas. 

a (A) 

c (A) 

9c = 90° — 2 tan ^ a/c (calc.) 

m 

None 

~ 6° 

« 8.0 

« 9.0 

6.7° 

0 

MgO 

5.2° 

8.11 

8.95 

5.63° 

m 

MgO 

5.1° 

8.14 

8.98 

5.62° 


TABLE II. Comparison of cubic-phase rotations 9 between experiment and prediction. The simple equation is accurate to 
within <1°. 














2 


The supercells were generated with ATAT’s MAPS code; then VASP fully relaxed the structure and calculated 
total energy. Varying supercell sizes require varying numbers of k-points; we use a P-centered mesh of at least 1000 
k-points per reciprocal atom, as implemented in ATAT. ATAT’s MAPS code then performs a cluster expansion to 
these energies. Then ATAT’s EMC2 code uses the same cluster expansion for phase diagram calculations. 

The structure generation, cluster determination, and effective cluster interaction (ECI) fitting was performed with 
the default settings of ATAT’s MAPS code. The calculation used a total of 192 ZnMnxGa 2 -x 04 supercells, of size 
< 42 atoms (« 460A ). The expansion’s cross-validation score is 1.1 meV/B-site (compare to energies of formation 
of 0-80 meV/B). A few additional structures relaxed to a high-energy peculiarity: 5 antiferromagnetic, 4 low-spin, 
2 with reoriented distortions (some or all octahedra distorted in y instead of z), and 2 anomalous cubic high-spin 
structures. All remained above the convex hull. When we restricted the antiferromagnetic and low-spin structures to 
a ferromagnetic high-spin configuration, the energy cost was only 4-12 meV/B. These are omitted from the expansion, 
because they simply add noise to the physics being fit in the expansion. Due to their high energy and minority, they 
cannot be expected to change the physics of the system. 

Eigure SI a) shows the quality of the fit for the various computed structures. Figure [sT|(b) shows the ECI as 
a function of cluster diameter. Intriguingly, the strongest interaction corresponds to the V 02 of Wojtowicz [ID], a 
pair repulsion between two adjacent atoms whose shared plane contains the JT-distorted direction. We suspect this 
interaction’s strength originates in the effect one cation’s <73 (JT-active) mode has on its neighbors’ octahedral modes 
and the associated energy cost. However, such a discussion is beyond the scope of this work. 
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FIG. SI. (a) Quality of fit in ATAT-generated structures. The CE energies and DFT energies match well, (b) Plot of ECI 
magnitudes as a function of distance. The picked CE used pairs, triplets, and quadruplets. 


We then computed the phase diagram using ATAT’s EMC2 code. We sampled temperatures from 100 to 1400 
K in steps of 100 K, with chemical potential (differences) /i from -30 to -1-30 meV/B in steps of 0.5 meV/B. For 
temperatures above 600 K, we additionally sampled y from -120 to -1-120 meV/B in steps of 0.5 meV/B and from 
-200 to -1-200 meV/B in steps of 5 meV/B. Results appear in Figurej^b) of the Letter. We verified the phase diagram 
with DFT-I-U, where Uun = 4 eV, using the same VASP -I- MAPS -I- EMC2 calculations. The phase diagram emerged 
qualitatively similar, although obviously the temperature scale differed. 

As mentioned in the Letter, long-range effects are notoriously difficult to capture with the CE. Later in the Letter, 


we analyze ZCO/ZMO slabs of thickness t layered in the direction k. Figure 
according to DFT and the above CE. 

From this figure, we conclude: 


S2 


compares the energies of these slabs 


1. The CE has a tendency to overmix relative to the exact result. This is because the phase-separation is due to 
the long-range JT effect; once the CE has difficulty capturing the full JT effect (see below), it will underestimate 
phase separation. Therefore, given that the CE phase diagram shows total phase separation, a fortiori this is 
the correct thermodynamic conclusion. 

2. The CE’s failure to precisely resolve these energies means the expansion is insufficiently robust to perform 
first-principles simulations of kinetics. A more accurate model will be necessary for that. 
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DFT / CE energies for multilayer slabs 



FIG. S2. Comparison of energies of multilayer ZGO/ZMO slabs, oriented in a particular direction, according to DFT and the 
cluster expansion (CE). 


The origin of the CE’s failure, as mentioned, lies in the long-range nature of the JT effect. This is corroborated 
by the long-range contact energies described in the Letter. Although the CE is a complete basis, in practice it is 
truncated after a finite number of terms. To successfully capture this, we would need to include many long-range and 
large (i.e. many-atom) clusters to the expansion. 

We attempted to capture these energetics with a variety of methods. 

First, we added these multilayer slabs to the CE’s list of known structures. However, the expansion still did not 
converge well. 

Second, we attempted including reciprocal-space clusters in our fit, via the mixed-basis cluster expansion [461152) . 
This method fits energies to clusters in reciprocal space that are defined by the structure factor: 

% = ^ CTj exp(ri • k) (SI) 

i 

where Oi = ±1 refers to the pseudospin of a particular site and Ti its location. However, this approach did not 
succeed either. The variety of necessary fc-space clusters, especially for the low-symmetry lii/amd crystal, makes 
this difficult. Additionally, the discrete phase boundaries make the A:-space expansion difficult. 

Third, we tried enhancing the fitting procedure with a compressive-sensing approach similar to that of Ref. 1531 and 
1541 This approach adds a penalty term for nonzero ECI in order to truncate the expansion after a few signihcant 
clusters. However, practically this penalty term takes the form of the ii norm (proportional to the ECI magnitude), 
rather than the £q norm (equal for all nonzero ECI). This leads to homogenization of ECI magnitudes, rather than 
physically expected decay of ECI with cluster size and distance. Despite attempts to compensate for this (similar to 
the reweighted normalization in [54] ). we could not find a convergent cluster expansion for this difficult system. 

Therefore, our CE is insufficiently robust for quantitative simulation, as evidenced by its failure to successfully 
predict the energies of ZGO/ZMO multilayer slabs. However, despite its tendency to overmix, it demonstrates the 
tendency of ZMGO to phase-separate even at high temperatures. 


SPIN-CROSSOVER CALCULATIONS 

It may be tempting to propose a spin crossover transition, either to explain the cubic structure at XMn < 0.25, or 
to cast aspersions on the GCA calculations. Perhaps Ga-doping raises the crystal held splitting (by shrinking crystal 
volume), moving the Mn’s d-band occupation from to Then the insignihcantly weak JT effect 

in the \t 2 g) manifold leads to an undistorted crystal strucutre. In fact, this spin-crossover transition was attributed 
to the tetragonal-to-cubic transition of ZMO at a pressure of 23 GPa [T71 HSl US]. Therefore, we further scrutinize 
GGA’s prediction of spin state and its reliability. 
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FIG. S3. ZMO’s spin-crossover transition. Energies (relative to the high-spin fully-relaxed ground state) are plotted as a 
function of fractional volume change. We show the checkerboard (CB) structure is in the high-spin regime according to GGA 
and a fortiori in experiment. 


However, it is clear that GGA is correct in predicting a high-spin state. As illustrated in Figure the reported 
nanocheckerboard lattice parameters place only a 5% (volumetric) strain on the Mn octahedra, whereas the ZMO 
phase transition occurs at a 10% (volumetric) strain [T7]; furthermore, GGA predicts the ZMO transition at only 
7% (volumetric) strain. It is obvious, then, that GGA has a tendency to overpredict the low-spin state, relative 
to experiment. Therefore, whereas GGA overpredicts the spin-crossover transition in ZMO, the absence of such 
a prediction in our structures—and the insufficient strain for such a transition—indicates that there is no spin- 
crossover transition in our structures or checkerboards. According to this, experiment should find high-spin Mn in 
ZnMno, 5 Gai. 504 , in agreement with GGA; no such experiments have been reported to our knowledge. 


SYMMETRY-ADAPTED HARMONICS 


In the Letter, we present formation energies for multilayer ZGO/ZMO slabs oriented in five directions. Here we 
detail the expansion to slabs oriented in an arbitrary direction. 

It is obvious from Figure |^b) that the slab formation energy of multilayer slabs of thickness t and cross-section A, 
oriented in direction k, can be expressed as £l(k, t) = A(ca(k) -|-te„(k)) , where Ca and Cy scale with area and volume 
(= At), respectively. The values of Ca,c„ for five directions are calculated from first principles as shown in Figure 

Sb). 

We now expand the coefficients c(k) as linear combinations of the symmetry-adapted harmonics, here the tetragonal 
(ZI 4 / 1 ) harmonics. These polynomials are orthonormal, where we take the inner product as: 


(/(k)l5(k)) 


f /(k)g(k) 




dip 


7r2 


fn{d,<p) 


cos 20 

cos 40 

e 

cos 4(^ sin'^ 0 

Order 



X* 

x^ 

Irrep. 

Alg 

Alg 


Big 

Ca (eV/A^) 

1 X 10"^ 

6 X 10“^ 

1 X 10"^ 

-9 X 10“'‘ 

c„ (eV/A") 

2 X 10"® 

-6 X10“^ 

5 X lO"'^ 

8 X 10“® 


(S2) 


TABLE III. Expansion of slab formation energies in symmetry-adapted harmonics. 

Table [III] shows the first four terms, along with a graphical representation, term order, and irreducible representation 
that the term transforms as. (The polynomial transforming as B^g also is of the order of x'^, but is a linear combination 
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of those listed.) The same table lists values for Cq and Cy, obtained by a least-squares fit of the four expansion terms 
to the five known energies. 


ENERGY DECOMPOSITION 

In the Letter, we describe the energy decomposition into strain and contact energy, and specifically its interpretation 
of the (Oil) preference. Here we present this decomposition in more detail. 

The formation energy of a ZnMnxGa 2 -x 04 structure (where XMn = x/2) is: 


^form — :^tot ^Mn^ZMO (1 ^Mn)^ZGO (^^) 

where all energies e are given per-B-site, and ezG(M)o refers to the energy of bulk ZG(M)0. The energy of formation 

is thus the difference between the energy and the “tie-line” connecting the ZGO / ZMO extrema. 

The energy of multilayer slabs are comprised of (a) bulk energies of the ZGO and ZMO slabs; (b) strain energy 
due to coherent lattice matching; (c) chemical binding energies; and (d) intracellular atomic displacements near the 
surface (k = 0 optical phonons). We combine these last two into a contact energy, noting that (unlike the strain) it 
must decay for infinitely thick slabs. We avoid the term “surface energy” because it can ambiguously refer to contact 
energy or total slab formation energy. 

Therefore, for two slabs of thickness t joined in direction k: 

^tot .^(^Mn^ZMO T (1 ^Mn)^ZGo) T T ^^contacti}^: (^^) 

We now seek to determine each of these terms from first-energy calculations. Consider the three structures shown 

in Figure!^ relaxed bulk (RB) of pure ZGO and ZMO, strained heterostructure (SH) of fully-relaxed adjacent slabs, 
and strained bulk (SB) of pure ZM(G)0 with lattice parameters set to those of SH. 



FIG. S4. Three types of ZMGO structures used to decompose slab formation energies. Arrows represent cross-sectional 
biaxial strain, calculated by fully relaxing the heterostructure. All calculations use periodic boundary conditions. See text for 
discussion. 

We can easily calculate the energy of each structure with DFT, using periodic boundary conditions. It is immediately 
apparent that 


Erb = N {xMnEzMO + (1 — a:Mn)l^ZGo) (S5) 

EsB = N{xMnEzMO + (1 “ a::Mn)l?ZGo) + NEstrain{^,t) (S6) 

ErH -^(:^Mn.^ZMO T (1 ^Mn).^ZGo) T (^; ^) “b (^^) 

It is then trivial to extract Egtrain and Econtact from DFT energies, as depicted in Figure [S4| 

Results are presented in the Letter (Fig. [^. For example, (001) layering contains little strain energy, because the 
JT-distorted [001] direction is perpendicular to the surface, so ZGO can achieve a cubic lattice and ZMO a tetragonal 
one. However, this is offset by significant contact energy, due to atomic rearrangements (interoctahedron coupling) 
near the surface. 

In the Letter we interpret the (Oil) preference with this decomposition. 
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